##Get objects
old_ws = ls()

####################
#Tables Intolerance#
####################

#intolerance data
intol_vars = c('live_in_village', 'live_in_neighbor', 'rents_room', 'marries_child', 'build_relg_estab')

#trust vars
trust_vars= c('trust_relg_more', 'trust_ethnicity_more', 'relg_candidate')

#safe vars
safe_vars = c('village_safe')


#Keep close elections, muslim respondents
useTable_ifls = mergedTable_ifls[(sample_ni_it) & (relg %in% 1)]
useTable_ifls[, w_i := 1 / .N 
                   , by = list(kab_code, dapil_number, election_cycle, cluster, IT_win)]

#Unique dapil-election cycle id
useTable_ifls[, dapil_election_id := paste(kabupaten, dapil_number, election_cycle, sep= ":")]

#Set relg_candidate NA for missing data
useTable_ifls[relg_candidate %in% c(6,9), relg_candidate := NA]

#observed?
useTable_ifls[, obs_1 := !(noone_present)]
useTable_ifls[, obs_2 := (spouse_present) | (adult_hh_present) | (adult_nh_present)]
useTable_ifls[, obs_3 := (adult_nh_present)]


############################
#Create Dependent Variables#
############################

#Recode for "yes" as 1, all other responses as 0. (Indicates yes to intolerant attitude)
useTable_ifls[, paste0(intol_vars, "_b") := lapply(.SD, function(x) x %in% c(1)), .SDcols = intol_vars]

#Reverse coding so higher values indicate intolerance.
useTable_ifls[, c(intol_vars) := lapply(.SD, function(x) (x * -1) + 5), .SDcols = intol_vars]

#Intolerance index (mean)
useTable_ifls[, intol_mean := rowMeans(.SD), .SDcols = paste0(intol_vars, "")]

#Intolerance index (mean)
useTable_ifls[, intol_mean_b := rowMeans(.SD), .SDcols = paste0(intol_vars, "_b")]

#All intolerant values
useTable_ifls[, intol_all := rowSums(.SD) == 5 , .SDcols = paste0(intol_vars, "_b")]

#Any intolerant values
useTable_ifls[, intol_any := rowSums(.SD) > 0 , .SDcols = paste0(intol_vars, "_b")]

#SUm intolerant positions
useTable_ifls[, intol_count := rowSums(.SD), .SDcols = paste0(intol_vars, "_b")]

#Recode for "yes" as 1, all other responses as 0. (Indicates yes to untrustworthy)
useTable_ifls[, paste0(trust_vars, "_b") := lapply(.SD, function(x) x %in% c(1)), .SDcols = trust_vars]

#Recode for so greater trst for ingroup is higher
useTable_ifls[, c(trust_vars) := lapply(.SD, function(x) x*-1), .SDcols = trust_vars]

#Trust index (mean)
useTable_ifls[, trust_mean := rowMeans(.SD), .SDcols = paste0(trust_vars, "")]

#All untrusting values
useTable_ifls[, trust_all := rowSums(.SD) == 3 , .SDcols = paste0(trust_vars, "_b")]

#Any untrusting values
useTable_ifls[, trust_any := rowSums(.SD) > 0 , .SDcols = paste0(trust_vars, "_b")]

#SUm untrusting positions
useTable_ifls[, trust_count := rowSums(.SD), .SDcols = paste0(trust_vars, "_b")]


##########################
#Intolerance and observed#
##########################

#Model specifications
lm_specs = expand.grid(x = c('IT_win*obs_3', 'IT_win*I(!obs_3)'), 
                       y = c(paste0(intol_vars, ""), paste0("intol", c('_mean', '_all', '_any'))),
                       stringsAsFactors = F
)


#Estimate
lm_models = Map(function(y,x) lm_loop(y,x, useTable_ifls, useTable_ifls[, w_i]),
                lm_specs$y, lm_specs$x)

#Standard errors
table_list = c(lm_models)
table_se = lapply(table_list, function(x) cluster_errors(x, useTable_ifls[, dapil_election_id]) %>% diag %>% sqrt)
table_cov = lapply(table_list, function(x) cluster_errors(x, useTable_ifls[, dapil_election_id]))
table_p = mapply(function(m,v) coeftest(m, vcov. = v)[4,4], table_list, table_cov, SIMPLIFY = F)

#Plot data
plot_data = data.table(observed = c("No", "Yes"), 
                       t = lapply(table_list, function(x) x$coef[2]) %>% unlist, 
                       se = lapply(table_se, function(x) x[2]) %>% unlist,
                       dv = c("Village",  "Neighbor", "Room", "Marriage" , "Building" ,"Mean", "All", "Any") %>% paste0("(",LETTERS[1:8],") ", .) %>% rep(each = 2),
                       p_diff = table_p %>% unlist)

plot_data[, dv := factor(dv, levels =c("Village",  "Neighbor", "Room", "Marriage" , "Building" ,"Mean", "All", "Any") %>% paste0("(",LETTERS[1:8],") ", .))]
plot_data[, lower := t - 1.96 * se]
plot_data[, upper := t + 1.96 * se]

text_data = plot_data[, list(dv, p_diff = paste0("p = ", sprintf("%01.3f", p_diff)))] %>% unique

pdf('./output/figures/figure_1.pdf', width = 11, height = 8.5)
p = ggplot(plot_data, aes(observed, t)) + geom_point() + 
  geom_errorbar(aes(ymin = lower, ymax = upper), width = 0) + facet_grid(. ~ dv) + 
  geom_hline(yintercept = 0, color = 'red') + 
  ylab("Estimated Effect of Islamist Win") + 
  xlab("Respondent Observed?") + theme_bw() +
  geom_text(
    data    = text_data,
    mapping = aes(x = -Inf, y = -Inf, label = p_diff),
    hjust   = -0.1,
    vjust   = -1
  )
print(p)
dev.off()

#Clean up
drop = setdiff(ls(), c(old_ws)) 
rm(list = drop)